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ABSTRACT 

We present methods and results from “21-cm Spectral Line Observations of Neutral Gas with the 
EVLA” (21-SPONGE), a large survey for Galactic neutral hydrogen (Hi) absorption with the Karl 
G. Jansky Very Large Array (VLA). With the upgraded capabilities of the VLA, we reach median 
root-mean-square (RMS) noise in optical depth of oy = 9 x 10 -4 per 0.42kms _1 channel for the 
31 sources presented here. Upon completion, 21-SPONGE will be the largest Hi absorption survey 
with this high sensitivity. We discuss the observations and data reduction strategies, as well as line 
fitting techniques. We prove that the VLA bandpass is stable enough to detect broad, shallow lines 
associated with warm Hi, and show that bandpass observations can be combined in time to reduce 
spectral noise. In combination with matching Hi emission profiles from the Arecibo Observatory 
(~ 3.5' angular resolution), we estimate excitation (or spin) temperatures (T s ) and column densities 
for Gaussian components fitted to sightlines along which we detect Hi absorption (30/31). We measure 
temperatures up to T s ~ 1500 K for individual lines, showing that we can probe the thermally unstable 
interstellar medium (ISM) directly. However, we detect fewer of these thermally unstable components 
than expected from previous observational studies. We probe a wide range in column density between 
~ 10 16 and > 10 21 cm -2 for individual Hi clouds. In addition, we reproduce the trend between cold 
gas fraction and average T s found by synthetic observations of a hydrodynamic ISM simulation by 
Kim et al. (2014). Finally, we investigate methods for estimating Hi T s and discuss their biases. 

Subject headings: ISM: clouds — ISM: structure — radio lines: ISM 


1. INTRODUCTION 

Star formation drives the evolution of galaxies, and 
the material from which stars form constantly cycles in 
and out of the interstellar medium (ISM). Within this 
process, neutral hydrogen (Hi) represents a crucial tran¬ 
sitionary stage between supernovae-expelled plasma and 
dense star-forming gas. Consequently, the presence and 
influence of radiative and dynamic physical processes in 
the ISM are imprinted in the mass distribution of Hi as 
a function of temperature and density. Observational 
constraints over the full range of temperature and den¬ 
sity are therefore necessary for realistic heating, cool¬ 
ing and feedback prescriptions within accurate models of 
star formation and galaxy evolution on all galactic scales 
(e.g. Bryan et al. 2007, Governato et al. 2010, Ostriker 
et al. 2010, Hill et al. 2012, Christensen et al. 2012). 

Analytical models of heating and cooling in the ISM 
predict the existence of two neutral Hi phases- the cold 
neutral medium (CNM) and the warm neutral medium 
(WNM)- each in thermal equilibrium (e.g., Pikel’ner 
1968, Field et al. 1969, McKee & Ostriker 1977). McKee 
& Ostriker (1977) and Wolfire et al. (2003) predict that 
the two phases should be mostly in thermal equilibrium, 
with the CNM dominating and the WNM comprising 
only a few percent of the total Hi column density. Over 
the range of pressures for which the two stable phases can 
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coexist, ~ 1000 < P/k < 8000 cm 3 K, we theoretically 
expect a volume density (n) and kinetic temperature (X*,) 
of (n,T K ) = (5-120 cm” 3 , 40-200 K) for the CNM, and 
re,TV) = (0.03-1.3 cm” 3 , 4100-8800 K) for the WNM 
Wolfire et al.l l2003lf . Within the high-density CNM, col¬ 
lisions between electrons, ions and other H atoms can 
thermalize the 21-cm transition so that T s = T U On the 
other hand, within the low-density WNM, collisions are 
not sufficient to thermalize the transition, and therefore 
it is widely expected that T s < (e.g. Field et al. 1958, 

Deguchi & Watson 1985, Liszt 2001). 

In addition, the presence of Hi in the intermediate, 
thermally unstable regime between the CNM and WNM 
has a range of theoretical implications. Dynamical pro¬ 
cesses are required to maintain this unstable material, as 
thermal pressure alone is not sufficient to force gas across 
extreme density boundaries (Cox et al. 2005). Shocks 
driven by supernovae into turbulent, magnetized gas cre¬ 
ate dense CNM clouds within a multiphase medium, 
where the fraction of thermally unstable gas is deter¬ 
mined by the supernova rate (Koyama & Inutsuka 2002, 
MacLow et al. 2005). Audit & Hennebelle (2005) showed 
that turbulent flow collisions produce unstable fractions 
of 10-30%, depending on the amount of turbulence in¬ 
jected. Modeling collisions of wind-blown superbubbles, 
Ntormousi et al. (2011) found most gas in cold clumps, 
with only 8-10% thermally-unstable WNM by mass. In 
another colliding flow model, Clark et al. (2012) showed 
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that the initial flow speed, including a time-dependent 
chemical network, has a huge effect on the resulting frac¬ 
tions of cold and warm gas. Clearly, observational con¬ 
straints of the WNM properties and the full temperature 
distribution over 10-10 4 K are important to constrain the 
role of self-gravity, turbulence, and large-scale gas flows 
in models of the ISM, as well as to guide initial conditions 
for numerical simulations. 

1.1. Previous Work 

Observationally, the hyperfine 21-cm transition of Hi 
provides an excellent tracer of the CNM and WNM. 
Given measurements of both Hi emission and Hi ab¬ 
sorption, it is possible to directly measure the excita¬ 
tion temperature (or spin temperature, T s ) and column 
density (N(Hl)) of CNM and WNM structures along 
the line of sight. Due to the high optical depth of the 
CNM, 21-cm absorption signatures are easy to detect, 
even with low sensitivity (e.g. Lazareff 1975, Dickey et 
al. 1977, Crovisier et al. 1978, Payne et al. 1978, Dickey 
et al. 1983, Braun & Walterbos 1992, Heiles & Troland 
2003, Kanekar et al. 2003, Mohan et al. 2004, Roy et al. 
2006, Begum et al. 2010). Extensive studies have been 
undertaken in the Milky Way and the CNM is easily 
detected all the way out to a Galactocentric radius of 
~ 25kpc (e.g., Dickey et al. 2009). 

However, the optical depth of the WNM is so low 
(t < 10~ 3 ) that extremely high sensitivity is required 
to detect it in absorption. Few studies have targeted 
individual detections of the WNM in absorption to di¬ 
rectly measure T s (e.g. Carilli et al. 1998, Dwarakanath 
et al. 2002, Murray et al. 2014), while many others es¬ 
timate WNM T s from upper limits (i.e. T/) or as line 
of sight averages in the presence of strongly absorbing 
CNM gas (Mebold et al. 1982, Heiles & Troland 2003, 
Kanekar et al. 2003, Kanekar et al2011, Roy et al. 2013a, 
Roy et al. 2013b). Therefore, the mass distribution of Hi 
in WNM-like conditions, i.e., at high T s and low n, is 
statistically unconstrained. 

Previous observational studies of CNM and WNM 
properties, although ubiquitous, typically do not reach 
sufficient sensitivity to detect the WNM or thermally 
unstable gas directly in Hi absorption. The Millennium 
Arecibo 21 Centimeter Absorption-Line Survey (Heiles 
& Troland 2003; HT03), comprised of 79 Hi absorp¬ 
tion and emission spectral pairs spread over the full 
Arecibo Observatory field of view, provided important 
constraints on the T s distribution of the CNM. A key re¬ 
sult from HT03 is that that 48% of the non-CNM spec¬ 
tral features had thermally unstable temperatures, with 
500 < T s < 5000 K. However, the high-T, measure¬ 
ments in the HT03 study are indirect (i.e., not based 
on direct detections of Hi absorption). Furthermore, 
single dish telescopes have the disadvantage that emis¬ 
sion fluctuations within the beam can contaminate ab¬ 
sorption measurements, while interferometers resolve out 
this large scale structure. In addition, the HT03 results 
could be biased by low sensitivity. Increased integration 
time on non-detection sightlines from HT03 reveal ab¬ 
sorption lines with r < 10 —3 , as shown by Stanimirovic 
& Heiles (2005). Previously-undetected, weak absorption 
lines with CNM-like T s will account for emission along 
the line of sight that was originally attributed to ther¬ 
mally unstable material, thereby reducing the estimated 


fraction of thermally unstable gas. Significant changes 
to the observed fraction of thermally unstable gas in the 
ISM will in turn have strong implications for the influ¬ 
ence of any radiative and dynamic processes responsible 
for creating and maintaining the instability. 

In a study of Hi absorption at the Karl G. Jansky Very 
Large Array (VLA), Begum et al. (2010) demonstrated 
that the VLA bandpass is stable enough to detect shal¬ 
low (peak r ~ 10 3 ), wide (full width at half maximum 
(FWHM) ~ 7 — 8kms _1 ) absorption lines. In combina¬ 
tion with Hi emission from the Arecibo Observatory, they 
identified individual absorption lines with spin tempera¬ 
tures in the thermally unstable regime, T s = 400 —900 K, 
in observations towards 12 bright background sources. 
However, only < 30% of their detected absorption lines 
had thermally unstable spin temperatures, which is much 
lower than the 48% unstable fraction reported by HT03. 
This emphasized the need for a larger interferometric 
study of Hi absorption at high sensitivity to further con¬ 
strain the fractions of gas in all Hi phases. 

Recently, Roy et al. (2013a,b) conducted an Hi ab¬ 
sorption survey of 35 continuum sources at the Wester- 
bork Synthesis Radio Telescope (WSRT), Giant Metre- 
wave Radio Telescope (GMRT) and Australia Telescope 
Compact Array (ATCA). Using the Leiden Argentine 
Bonn (LAB; Kalberla et al. 2005) survey for Hi emission, 
they found at least 28% of the absorption-detected Hi in 
the thermally unstable regime (Roy et al. 2013b). How¬ 
ever, this value is based on only 13 of their sources and 
the conservative assumptions that all detected CNM has 
T s = 200 K and all non-detected WNM has T s = 5000 K, 
and not on direct measurements of the absorbing and 
emitting properties of the unstable gas. 

1.2. 21-SPONGE 

To measure the physical properties of all neutral Hi in 
the ISM, we are conducting a large statistical survey, “21- 
cm Spectral Line Observations of Neutral Gas with the 
EVLA” (21-SPONGE), to obtain high-sensitivity Milky 
Way Hi absorption spectra using the VLA. The recently 
upgraded capabilities of the VLA allow us to achieve me¬ 
dian RMS noise levels in optical depth of o> = 9 x 10~ 4 
per 0.42 kms' 1 channel, which are among the most sen¬ 
sitive observations of Hi absorption to date. Currently 
31 sources are complete after over 200 hours of observing 
time. We have a very high detection rate so far, and we 
detect Hi absorption in the direction of 30/31 sightlines. 

In Murray et al. (2014) we applied a spectral stacking 
technique to the initial one third of the 21-SPONGE sam¬ 
ple (19 Hi absorption lines) and detected a widespread 
WNM component at high significance. The spin tem¬ 
perature of this component, equal to 7200+ 4 ®qo K (68% 
confidence), is higher than predictions based on colli- 
sional excitation alone (< 4400K, Liszt 2001), implying 
that additional excitation mechanisms must be prevalent 
in the ISM. Lya scattering provides a likely source for 
this additional excitation. However, the required den¬ 
sity of Lya photons necessary to explain our result is 
significantly higher than expectations from recent theo¬ 
retical and numerical studies (e.g. Liszt 2001, Kim et 
al. 2014). For example, in their Galactic simulations, 
Kim et al. (2014) used a fixed Lya photon number den¬ 
sity of 10 -6 cm -3 and found T s = 4000 — 5000 K for the 
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TABLE 1 

Survey Comparison 


Survey 

Galactic 
latitude (°) 

(i) 

Angular res. 
(arcmin) 
(2) 

Velocity res. 
(kms -1 ) 
(3) 

Sensitivity 

(4) 

Number 
of spectra 
(5) 

Telescopes 

(absorption, emission) 

(6) 

Absorption targeted: 


21-SPONGE 

|*| > 3.7 

4 

0.5 

0.0006 

58 

VLA, Arecibo 

Roy et al. 2013 

|b| > 1.6 

36 

1.0 

0.0005 

35 

WSRT/GMRT/ATCA, LAB 

Heiles & Troland 2003 

all |b| 

4 

0.16 

0.002 

78 

Arecibo 

Stanimirovic et al. 2014 

—32 <b <-8 

4 

0.16 

0.002 

27 

Arecibo 

Mohan et al. 2004 

|6| >15 

36 

3.3 

0.005 

102 

GMRT, LDS 

Emission mapped: 


VGPS b 

|6| <1.3 

1 

1.56 

0.025-0.125 

113 

VLA, GBT 

CGPS C 

-3.6< b <15 

1 

1.32 

0.023-0.115 

364 

DRAO, DRAO 26m 

SGPS d 

|6| <1.5 

2 

1.0 

0.02-0.1 

96 

ATCA, Parkes 

GASKAP 6 

|6| <10 

0.5 

1.0 

0.02 

>1000 

ASKAP 


a : Median RMS noise in Hi optical depth per lkms 1 channel. 

b : IStil et ahl (120013 ) ; c : ITavlor et al.l (120031 1 ; d : IMcClure-Griffiths et all 320051 1 ; e : IDickev et al.l 320131 1 


WNM. Understanding the excitation processes of Hi, as 
well as the strength and topology of the Lyct radiation 
field, clearly require future studies. 

In Table [1] we compile properties of other single dish 
and interferometric surveys for Hi absorption for com¬ 
parison, including: (1) estimates of the area covered in 
Galactic latitude (°), (2) the angular resolution in Hi 
emission (arcmin), (3) the velocity resolution (kms -1 ), 
(4) the Hi optical depth sensitivity (oy) per lkms -1 
channels, (5) the number of Hi absorption spectra, and 
(6) the telescopes used for Hi absorption and emission. 
Those studies which targeted particular sources for mea¬ 
suring Hi absorption at a desired sensitivity are classi¬ 
fied as “Absorption targeted”, and those studies which 
mapped large spatial areas in Hi emission and then ex¬ 
tracted absorption from continuum sources within the 
map are classified as “Emission mapped”. We note that 
the 4 surveys in the latter category were all conducted 
at low-|6| and therefore their spectra are more complex 
than those at high-latitude in the targeted absorption 
line surveys. 

21-SPONGE is one of the most sensitive and extensive 
absorption line surveys ever undertaken. We emphasize 
that our analysis will improve on previous work at sim¬ 
ilar sensitivity, as we estimate physical properties of in¬ 
terstellar Hi using emission data obtained with spatial 
resolution closer to the ~ 1" interferometric absorption 
resolution from the VLA than previous studies. For ex¬ 
ample, Roy et al. (2013) use emission from the ~ 36'-wide 
LAB beam, and we use emission data from the Arecibo 
Observatory, whose beam area (~ 3.5' x 3.5') is smaller 
than the LAB beam by two orders of magnitude. In addi¬ 
tion, 21-SPONGE will occupy a unique space in terms of 
sensitivity and size for a long time, given that future Hi 
absorption surveys with the Australian Square Kilome¬ 
ter Array Pathfinder (ASKAP) and the Square Kilometer 
Array (SKA) will not likely seek such high sensitivity in 
optical depth (e.g. McClure-Griffiths et al. 2015). 

In this paper, we present the current progress of 21- 
SPONGE. In Section 2, we describe the observations and 
data reduction; in Sections 3 we discuss our methods of 
spectral line analysis via Gaussian fitting; in Secion 4 we 
show example decompositions for select sources; in Sec¬ 
tion 5 we calculate temperatures and column densities 
for all individual component fits; in Section 6 we com¬ 


pare our results with previous observational surveys; in 
Section 7 we compare our results with synthetic obser¬ 
vations of a hydrodynamic Galaxy simulation by Kim et 
al. (2014); in Section 8 we discuss the pros, cons and bi¬ 
ases of various methods for calculating T s from Hi emis¬ 
sion and absorption; and in Section 9, we summarize the 
results. 

2. DATA PROCESSING 
2.1. Observations 

VLA observations for 21-SPONGE began in February 
2011, and targeted 58 sources from the NR AO/VLA Sky 
Survey (NVSS) catalog with flux densities at 1.4 GHz 
(<Si. 4 GHz) > 3 Jy (Condon et al. 1998). We chose most 
sources to lie at high galactic latitude (|&| > 10°) to 
avoid complicated Hi profiles associated with the Galac¬ 
tic plane, and also to have angular sizes less than 1 ' to 
avoid resolving substantial continuum flux. The range of 
source parameters allows us to conduct observations in 
all array configurations and configuration-moves except 
for D-array (most compact configuration). Currently, 53 
sources have received full integration time and 31 of these 
53 have been fully reduced and analyzed. Table [2] con¬ 
tains information about the completed observations. 

All VLA observations use simultaneously three sep¬ 
arate, standard L-band configurations, each with one 
dual-polarization IF of bandwidth 500 kHz and 1.95 kHz 
per channel spacing. One IF is centered at 1.420408 GHz 
(Hi line rest frequency, “standard”), one at 1.421908 GHz 
(1.5 MHz, or about 300 km s -1 , higher than the Hi 
rest frequency, called “high”) and one at 1.418908GHz 
(1.5MHz lower, called “low”). We use the high and low 
configurations for observing our bandpass calibrators to 
avoid contamination by local Hi in the direction of all 
calibrators at the Hi rest frequency. Because we nor¬ 
malize our solution with respect to the continuum level, 
the absolute phase change associated with the frequency 
switching method is not an issue. We use the “standard” 
configuration for observing target sources and measur¬ 
ing Hi absorption, allowing for a velocity coverage of 
107.5 kms -1 with 0.42 km s -1 channels, corresponding 
to a velocity resolution of 0.5kms -1 (Rohlfs & Wilson 
2004). Relative flux calibration is performed via self¬ 
calibration on the target source and phase and ampli¬ 
tude calibration are performed by observing a nearby 
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TABLE 2 

VLA Observation Information 


Source 

(name) 

RA(J2000) 

(hh:mm:ss) 

Dec (J2000) 
(dd:mm:ss) 

1 

(°) 

b 

(°) 

Si. 4 GHz 

(Jy)° 

a T 

(10“ 3 ) 6 

Synthesized Beam 
(arcsec 2 ) 

Time 

(hrs) c 

4C32.44 

13:26:16 

31:54:10 

67.240 

81.049 

4.6 

1.5 

2.8x1.3 

6.6 

4C25.43 

13:30:37 

25:09:11 

22.464 

80.991 

6.9 

0.9 

2.7x1.2 

2.9 

3C286 

13:31:08 

30:30:33 

56.526 

80.676 

14.9 

0.7 

4.1x3.0 

2.8 

4C12.50 

13:47:33 

12:17:24 

347.220 

70.173 

5 

1.3 

4.2x1.5 

2.8 

3C273 

12:29:06 

02:03:05 

289.945 

64.358 

54.9 

0.6 

7.8x4.2 

1 

3C298 

14:19:08 

06:28:35 

352.159 

60.667 

6 

0.8 

2.3x1.4 

2.8 

4C04.51 

15:21:14 

04:30:19 

7.290 

47.748 

4 

2.9 

12.2x4.7 

5.6 

3C237 

10:08:00 

07:30:16 

232.117 

46.627 

6.5 

0.9 

6.5x4.4 

5 

3C225A 

09:42:15 

13:45:51 

219.866 

44.025 

4.4 

1.5 

3.9x1.6 

6 

3C225B 

09:42:15 

13:45:49 

220.010 

44.007 

4.4 

3.0 

3.9x1.6 

6 

3C345 

16:42:59 

39:48:37 

63.455 

40.948 

7 

1.3 

3.4x1.4 

2.5 

3C327.1 

16:04:45 

01:17:51 

12.181 

37.006 

4.1 

1.4 

6.1x4.8 

4.7 

3C147 

05:42:36 

49:51:07 

161.686 

10.298 

22.9 

0.4 

4.4x3.8 

1.1 

4C33.48 

19:24:17 

33:29:29 

66.388 

8.373 

3.8 

2.8 

4.8x1.8 

5.7 

3C154 

06:13:50 

26:04:36 

185.594 

4.006 

5.2 

0.9 

14x13 

3.5 

3C410 

20:20:06 

29:42:12 

69.210 

-3.768 

10 

1.9 

9.9x5.5 

2.5 

B2050+36 

20:52:52 

36:35:35 

78.858 

-5.124 

5 

2.3 

4.3x2.0 

5 

3C409 

20:14:27 

23:34:52 

63.397 

-6.120 

14 

3.3 

1.6x1.4 

3.3 

PKS0531+19 

05:34:44 

19:27:21 

186.761 

-7.109 

7 

0.5 

1.3x1.1 

1.9 

3C111 

04:18:21 

38:01:35 

161.675 

-8.821 

4.3 

0.7 

23x23 

5.1 

3C133 

05:02:58 

25:16:24 

177.725 

-9.914 

5.7 

1.8 

1.3x1.1 

5.1 

3C138 

05:21:09 

16:38:22 

187.403 

-11.346 

9 

0.9 

15x5.2 

2.2 

3C123 

04:37:04 

29:40:13 

170.581 

-11.662 

47 

0.7 

19x5.5 

2.2 

3C433 

21:23:44 

25:04:10 

74.475 

-17.693 

12 

2.4 

9.7x5.7 

3 

3C120 

04:33:11 

05:21:15 

190.373 

-27.397 

3.4 

1.1 

4.7x4.3 

7 

3C48 

01:37:41 

33:09:35 

133.961 

-28.720 

15.7 

0.7 

1.3x1.2 

2.8 

4C16.09 

03:18:57 

16:28:32 

166.633 

-33.598 

8 

0.8 

1.3x1.2 

1.5 

3C454.3 

22:53:57 

16:08:53 

86.108 

-38.182 

11 

0.9 

7.5x3.9 

2.5 

J2232 

22:32:36 

11:43:50 

77.436 

-38.582 

7 

0.7 

5.2x4.3 

2.1 

3C78 

03:08:26 

04:06:39 

174.857 

-44.514 

5.7 

4.1 

4.6x3.7 

5.2 

3C459 

23:16:35 

04:05:17 

83.038 

-51.285 

4 

0.8 

5.7x4.5 

5 


“: Condon et al. (1998). 

b \ RMS noise in r per 0.42 kins -1 channel, measured off-line. 
c : Total on-source time, not including calibration overheads. 

VLA calibrator source. 

For bandpass calibration, we observe the strong cal¬ 
ibrators 3C286, 3C48 or 3C147 in both high and low 
instrument configurations following each on-source scan. 
For each full-length observing session (between 3 and 5 
hours) we allocate up to 80% overhead for the high and 
low bandpass observations. However, if the source was 
stronger than ~ 7 Jy, it was observed as its own band¬ 
pass calibrator to conserve slewing time over the course of 
the observation, which is particularly beneficial for filler- 
project length observing files (maximum 1 hour total). 

2.2. Data reduction 

All data sets were reduced using the Astronomical Im¬ 
age Processing Systeir0 (AIPS). We decided to use AIPS 
instead of CASA0 to take advantage of the suite of band¬ 
pass analysis tasks unique to AIPS when we began the 
survey. During observations conducted in configuration 
moves involving D-array, all baselines shorter than 300 m 
were excluded to avoid contribution from partially re¬ 
solved Hi emission. After initial flagging of all observa¬ 
tions, the task BPLOT was used to examine initial band¬ 
pass solutions for the high and low bandpass calibration 
observations in detail for each antenna. These solutions 

2 http://www.aips.nrao.edu 

3 http://casa.nrao.edu/ 


were then combined to produce a final bandpass solu¬ 
tion. For all sources, we consistently reach noise levels 
in normalized bandpass amplitude profiles of < 10 -3 per 
0.42kms _1 channels. 

In addition to a linear slope in the bandpass solution, 
we found a stable, periodic ripple which was revealed 
only because of the extreme sensitivity of our observa¬ 
tions. The ripple is caused by the finite impulse response 
(FIR) filters applied to the data prior to correlation, and 
is stable in amplitude, phase and time. Appendix A con¬ 
tains information on the properties and stability of the 
bandpass structure. 

Furthermore, we demonstrated that the VLA band¬ 
pass is stable enough to combine bandpass observations 
conducted at different times or on different days. This is 
especially helpful for filler observations (max total obser¬ 
vation time <1 hour). We are able to combine bandpass 
observations between separate days, thereby increasing 
the total bandpass observation time. The noise level in 
the bandpass solution decreases with added time by in¬ 
dex of —0.44 ± 0.02, which is nearly the theoretically 
expected index of —0.5 (i.e. by l/y/tz/ti for combined 
integration time t 2 and initial integration time t{). Ap¬ 
pendix A contains further details on the noise improve¬ 
ment due to time averaging in the cases of additional 
datasets. 

Following amplitude and phase calibration on the tar- 
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get dataset, we apply the bandpass solution from the 
combined bandpass calibration observation. We then 
perform self-calibration on the target dataset by isolating 
the target source continuum, constructing an image, and 
using it to calibrate the continuum dataset. After apply¬ 
ing the self-calibration solution, we repeat the process 
until the signal-to-noise level in the resulting continuum 
image no longer improves significantly. For most sources, 
this occurs after only two iterations. Returning to the 
original target dataset, the continuum contribution to 
each source is determined by fitting a linear model to 
line-free channels. We subtract this from the source vis¬ 
ibilities in all channels using the AIPS task UVLSF. We 
then use the task CVEL to correct source velocities for 
Earth’s rotation and motion towards the Local Standard 
of Rest (LSR) within the Solar System^! 

The final calibrated data cubes are constructed using 
the AIPS task IMAGR, and are cleaned to three times 
the background noise in a test image of one channel. 
For each cube, the pixel sizes are calculated to be 4 
times smaller than the beam size to properly sample the 
beam. Finally, we extract the absorption spectrum from 
the central pixel of each image, as in most cases, the 
sources are unresolved by the observations. However, 
of the 31 completed sources, 5 were resolved (3C327.1, 
3C154, 3C111, 3C133, 3C123). To extract the spectra 
from these 5 data cubes, we average all spectra over the 
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Fig. 1.— Comparing integrated optical depth values between 
21-SPONGE and Roy et al. 2013 (R13; black stars), 21-SPONGE 
and HT03 (green diamonds) and 21-SPONGE and Stanimirovic Sz 
Heiles (2005) (SH05; filled green diamond). The dotted line shows 
y=x to highlight the agreement between the majority of points. 
Sources which strongly deviate from y=x are indicated by their 
source names. Errors on the measurements are included, but are 
smaller than the symbol sizes unless visible. 


4 We initially applied CVEL before continuum subtraction 
(UVLSF). However, we observed considerable Gibbs ringing in the 
edge channels of our reduced spectra. This forced us to Hanning 
smooth all datasets, reducing our velocity resolution by a factor of 
2. After investigating the origin of the Gibbs ringing, we found that 
the ringing phenomenon arose in the order of CVEL and continuum 
subtraction, and reversing their order of application eliminated the 
ringing effect. 


extent of the source. 

The achieved sensitivities for the final data products 
are listed in Table [2l and the quoted values for oy are 
per 0.42 km s -1 channel. The median RMS noise values 
are oy = 9 x 10 -4 (mean 1 x 10 -3 ) per 0.42kms _1 chan¬ 
nel, and oy = 6 x 10~ 4 (mean 8 x 10“ 4 ) per lkrns^ 1 
channel. Several sources have high noise due to the fact 
that they were resolved by the VLA, and we plan to in¬ 
clude additional data in these directions to reduce the 
noise for the final data release. 

2.3. Optical Depth Comparison 

To verify the accuracy of our absorption spectra, 
we compare integrated optical depths with two previ¬ 
ously published surveys: the Arecibo Millennium Survey 
(HT03) composed of 79 absorption spectra with match¬ 
ing expected emission spectra from the Arecibo Obser¬ 
vatory, and 35 high-sensitivity absorption spectra from 
the GMRT, WSRT and ATCA by Roy et al. (2013). We 
compare the sources overlapping between 21-SPONGE 
and HT03 (green diamonds) and between 21-SPONGE 
and Roy et al. (2013) (black stars) in Figure Q] 

The three catalogs generally agree well in terms of inte¬ 
grated optical depth. In Figure [TJ we indicate the source 
names of those which differ most between the catalogs. 
In the cases of 4C25.43 and 3C345, and Roy et al. (2013) 
have better sensitivity (ay = 5 x 10 -4 and 5 x 10~ 4 
whereas we achieve oy = 13 x 10~ 4 and 9 x 10 -4 re¬ 
spectively), so our f r(v)dv values have larger error bars. 
For source 3C273, the HT03 baseline is not perfectly flat, 
with the consequence that their f r(v)dv value is higher 
than ours. HT03 did not measure any Hi absorption 
above the la limit towards 3C286, so we compare our 
optical depth spectrum to the higher-sensitivity result 
from Stanimirovic & Heiles (2005; SH05) (filled green di¬ 
amond), and observe that our integrated values are con¬ 
sistent. In addition, we agree very well with the Roy 
et al. (2013) value for 3C286. Towards source 4C32.44, 
HT03 have lower sensitivity and so their value has larger 
error bars. Finally, the HT03 profile towards 3C225A, 
a resolved double-lobed radio galaxy, contains “positive 
absorption” due to an artifact caused by the Hi emis¬ 
sion fluctuations and large uncertainties in estimating 
Hi emission profiles. 

Remaining residual differences that exist could possi¬ 
bly be caused by fluctuations in absorption profiles due 
to transient clouds, as has been well-studied in the cases 
of 3C147 (Lazio et al. 2009) and 3C138 (Brogan et al. 
2005, Roy et al. 2012). These studies find differences in 
optical depth between observational epochs (spaced by 
several years) of up to ~ 0.5. 

Overall, the excellent agreement between 21-SPONGE 
absorption and HT03, Roy et al. (2013) and Stanimirovic 
& Heiles (2005) absorption demonstrates that our ob¬ 
serving and reduction strategies are sound. It is also 
encouraging to see that in the majority of cases, 21- 
SPONGE observations agree well with HT03 and Sta¬ 
nimirovic & Heiles (2005). This proves that Arecibo - 
a single dish telescope - can reliably produce accurate 
Hi absorption profiles in the directions of strong sources 
(Si. 4 GHz > 3 Jy). 

2.4. Matching Emission Profiles 
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We obtained matching Hi emission profiles along most 
sightlines from the 305-m Arecibo Observatory as part 
of project a2770. For the observations, we used the L- 
wide receiver to simultaneously observe Hi at 1420 MHz 
and three OH lines (1665, 1667 and 1720 MHz) with 
two linear polarizations. We achieve a velocity reso¬ 
lution of 0.16 kms -1 over 2048 channels between —164 
and 164kms _1 for the Hi observations. The ~3.5’ angu¬ 
lar resolution of the Arecibo telescope at these frequen¬ 
cies complements the VLA observations by minimizing 
the effects of mismatched beam sizes on interpreting the 
absorption and emission spectra. It would require pro¬ 
hibitively long integration times to acquire emission pro¬ 
files at similar sensitivity from the VLA in D array. 

Following the methods of HT03, we constructed an 
“expected” Hi emission profile (T exp (v)) toward each 
source by interpolation using a pattern of 16 observed 
off-source positions. The expected profiles represent the 
profile one would observe if the continuum source were 
turned off. The profile construction via least squares fit 
is described fully by HT03, and was also implemented 
by Stanimirovic & Heiles (2005) and Stanimirovic et 
al. (2014). The error in T exp (v) as a function of veloc¬ 
ity, cr (u), is also computed in this process based on 
the difference in system temperature and spatial offset 
between on-source and off-source pointings (HT03, Sec¬ 
tion 2.7). As discussed by Stanimirovic et al. (2014), we 
used a simpler method than HT03 and excluded fine 
tuning for gain variations under the assumption that 
correctly accounting for them requires accurate knowl¬ 
edge of the spatially-varying gain and beam shape. In 
addition, we used a second-order Taylor expansion to 
construct T exp (v), which, as found by Stanimirovic et 
al. (2014), is a noisier but more accurate approach than 
using a first-order Taylor expansion. We then divided 
the expected profiles by a total beam efficiency factor of 
0.85 (P. Perillat private communication) to convert the 
expected profiles to brightness temperature. 

To compare our results, we scaled the expected emis¬ 
sion profiles from HT03 to our new Arecibo expected 
emission profiles by their velocity integrals and found 
that applying a median beam efficiency factor of 0.81 
to the HT03 profiles produced the best agreement with 
our beam efficiency-scaled profiles from a2770. For each 
source, we then selected the scaled emission profile with 
the best sensitivity. For 11/31 sources we use the HT03 
expected emission profiles, and for 16/31 sources we use 
a2770 data. The remaining 4 sources do not have Arecibo 
emission profiles because they are on the edge of the 
Arecibo field of view, and so we use emission spectra 
from the Leiden Argentine Bonn (LAB) survey for these 
sources (Kalberla et al. 2005). 

2.5. Stray Radiation 

Due to the complex beam pattern of the Arecibo tele¬ 
scope, our emission profiles are likely cotaminated at 
some level by radiation from higher order side-lobes, 
known as stray radiation. To estimate the level of this 
contamination, we use additional data from the Galactic 
Arecibo L-band Feed Array Survey in Hi (GALFA-Hi; 
Stanimirovic et al. 2006, Peek et al. 2011). GALFA-Hi 
is an all-sky survey for Hi emission in the Galaxy using 
the seven beam-array ALFA receiver. Peek et al. (2011) 
smoothed GALFA-Hi data cubes to the ~ 36' angular 


resolution of the stray radiation-corrected LAB survey 
(Kalberla et al. 2005) and found that the differences gen¬ 
erally fall within typical ler errors (see Figure 12 in Peek 
et al. 2011). Stray radiation signatures appear as broad, 
wing-like features at high velocities, which are easy to 
isolate visually from the narrower, peak-like differences 
due to mismatched beam shapes between the LAB sur¬ 
vey and the smoothed GALFA-Hi data cubes. 

To compare our results for all sources covered by 
GALFA-Hi to date, we constructed expected GALFA- 
Hi emission profiles around each source by averaging 
5x5 pixel 2 regions around each source. The central 
pixels of this average are likely affected by absorption, so 
we removed the central 3x3 pixel 2 region which corre¬ 
sponds to roughly one beam area. We then subtracted 
our T exp (v) profiles from the expected GALFA-Hi pro¬ 
files and found that the differences fall mostly within 3 ct 
errors in T exp (v ), which is very encouraging. The best 
method of removing the stray radiation contamination 
in lieu of modeling the complex, zenith-angle-dependent 
beam pattern of the Arecibo telescope would be to boot¬ 
strap the GALFA-Hi data to the LAB survey and simply 
remove the differences from our profiles (e.g. as discussed 
by Peek et al. 2011). However, applying this method 
necessarily injects noise from the LAB survey into the 
higher-sensitivity GALFA-Hi profiles. In addition, the 
method assumes that the contamination to GALFA-Hi 
(observed with ALFA) is the same as to our T exp (v ) pro¬ 
files (observed with L-wide), which may not necessarily 
be the case. In the future we plan to investigate alter¬ 
native methods for modeling the stray radiation contri¬ 
bution to our emission profiles without degrading our 
sensitivity unnecessarily. 

3. ANALYSIS 
3.1. Gaussian Fitting 

To estimate physical properties of individual interstel¬ 
lar clouds along the line of sight, we decompose the Hi 
spectra into Gaussian functions (e.g. HT03). Although 
the CNM dominates the absorption spectra, by achiev¬ 
ing extremely high sensitivities we aim to detect shallow, 
wide signatures of warm gas directly in absorption. 

Based on the two-ph ase ISM model, discussed in detail 
bv lDickev et all (1200311 and HT03, both CNM and WNM 
contribute to the expected profile, T exp (v). We denote 
the CNM and WNM components detected in absorption 
by the subscript “ABS” and the WNM components de¬ 
tected in emission by the subscript “EM”. Thus T exp (v) 
is given by, 

Texp (v) = T b ,abs(v) + T b ,em(v), (1) 

where T B , ABs(f) is the brightness temperature due to 
absorption-detected gas, and Tb.em('d) is the brightness 
temperature due to non-absorbing gas detected only in 
emission. 

For each absorption spectrum, we identify a number 
(N) of Gaussian components using a least-squares fit (de¬ 
scribed in Section 3.3), so that the absorption spectrum 
t(v) is given by, 

iV-l 

t(v) = T 0 , n e- [ ( v - Vo ^ /SVn] \ 

n =0 


(2) 
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where (to, vo, n , Sv n ) are the (peak optical depth, central 
velocity, 1/e width) of the n th component in the fit. 

Following the fit to the absorption spectrum, and as¬ 
suming each fit component is independent and isothermal 
with spin temperature T Sj71 , the brightness temperature 
contribution from the N Gaussian functions detected in 
absorption to the expected emission profile is given by, 

N-l 

T B , abs(v) = T ^ 1 - (3) 

n =0 

Here, the subscript m and the associated optical depth 
T m (v) run through the M number of CNM clouds lying 
in front of cloud n. 

We model the remaining contribution to the expected 
profile by WNM not detected in absorption using a set 
of K independent Gaussian functions. For each k “ com¬ 
ponent, we assume that a fraction F\ lies in front of all 
N absorption components, and thus a fraction (1 — F k ) 
is absorbed by the intervening gas. The emission bright¬ 
ness temperature from these K non-absorbing Gaussian 
components is given by, 


K -1 

T b , em(v) = ^[F fc +(l-F fc )e- T W]xTo, fe e- [ ^-^ fe )/^l 2 , 
k =0 

(4) 

where (To.fe,are the (peak in units of bright¬ 
ness temperature, central velocity, 1/e width) of the k th 
emission component. For a given set of absorption com¬ 
ponents, the expected profile is fit for the spin tempera¬ 
tures of the N absorption components and the Gaussian 
parameters of the K emission components. This fit is 
repeated for all permutations of absorption components 
along the line of sight, for every value of Fk for each 
emission component. Following the method of HT03, we 
allow Fk to be 0, 0.5, or 1. The values of Fk affect the 
derivation of spin temperatures, but any finer tuning of 
the values produces results that are difficult to distin¬ 
guish statistically. Overall, there are N\ x 3 A permuta¬ 
tions for each line of sight. The final fit selected is the 
one with the smallest residuals, and the final spin tem¬ 
peratures are calculated by a weighted average over all 
trials, where the weights are equal to the inverse of the 
variance of the residuals to the T exp (y) fit (HT03). The 
error in the final spin temperature is estimated from the 
variation in T s with Fk, as described in Section 3.5 of 
HT03. 

The expected emission profile (Equation [TJ has been 
baseline-corrected, meaning that the contribution from 
diffuse radio continuum emission, including the cosmic 
microwave background (CMB) and Galactic synchrotron 
emission, has been removed. This value, T s k y , is equal to 
2.725 K from the CMB plus an estimate of the Galactic 
synchrotron background at the source position. Given 
that the spectral index of the synchrotron background is 
about 2.7, we divide the 408 MHz brightness tempera¬ 
ture taken from the Haslam et al. (1982) Galactic survey 
by a factor of (1420/408) 2 '' to estimate the synchrotron 
contribution at 1420 MHz. The 21-SPONGE sources are 
mostly at high latitude, and therefore the values of T s k y 
are generally small, ranging between 2.76 K and 2.85 K. 
To account for its presence, we add T s k y back to T exp (v) 


before performing the fit. 

3.2. Absorption Noise Spectra 

The noise in Hi absorption spectra is not constant as 
a function of velocity. Hi emission at Galactic velocities 
will increase the system temperature of the VLA anten¬ 
nas, producing different noise properties across the ab¬ 
sorption spectra. To estimate the “noise spectra”, a T (v), 
of the VLA absorption profiles, we follow the methods 
discussed by Roy et al. (2013a). The spectral noise is 
composed of (1) on-source noise, a on (v) (varies with ve¬ 
locity) and (2) off-line noise from the frequency-switched 
bandpass solution, ctbp (constant with velocity). We first 
estimate the RMS noise in off-line channels, ay (see Ta¬ 
ble 2), which includes contributions from both a on (v) 
and ctbp- We then scale ay by ((Tb(u) + T sys )/T sys ) to 
calculate a on (v). The antenna used to produce the LAB 
survey is of similar size to an individual VLA antenna, 
so we use Hi emission from the LAB survey to estimate 
the brightness temperature as a function of frequency 
for each source (Tb(v)) and we assume that the system 
temperature of the VLA in the absence of Hi emission 
is T sys ~ 25 K. Finally, we calculate the noise spectrum 
by: oy(u) 2 = a on (v ) 2 + a PP (also see discussion by Roy 
et al. 2013a). Examples of these noise spectra are dis¬ 
played as dotted lines above and below the residuals of 
the Gaussian fits in Figures [H [3] and [4] We use the 
noise spectra to visually inspect the fit residuals, not to 
quantify the goodness-of-fit. 

3.3. Component Selection: Absorption 

The process of Gaussian fitting has many inherent un¬ 
certainties, including the determination of number of fit 
components. The noise level in the fit residuals can be 
used as a goodness-of-fit indicator, but it cannot eas¬ 
ily distinguish between different numbers of components, 
because adding more components will always reduce the 
residual noise level without necessarily bringing the fit 
any closer to an accurate depiction of reality. 

To determine the best-fit number of Gaussian com¬ 
ponents to the absorption spectra, we begin by using a 
reduced chi-squared (x^ed) test- We calculate Xr e d using: 


2 _ 1_ (0(v)-M(v)) 2 

Wed v 2^ 


o-f 


(5) 


where 0(v) are the observed data, M(y) are the fitted 
model data, ay is off-line RMS noise in the absorption 
spectrum (from Table 2) and v is the number of degrees 
of freedom (DOF: equal to number of measurements (256 
spectral channels) minus the number of fitted parame¬ 
ters, which includes the height, width and center of each 
Gaussian function and a constant continuum offset). 

For each profile, we would like to include the minimum 
number of components that brings \ 2 ed near to 1.0. How¬ 
ever, in many cases, the difference in x^ed between fits 
of different numbers of components is too small to confi¬ 
dently determine the best fit. As an example, we present 
the case of 3C123 in Figured The central component of 
the absorption spectrum is very strong, and is poorly fit 
by a single Gaussian component. Including an additional 
broad component at the base of the central component 
does improve the fit, and including additional narrow 
lines improves the fit further. 
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Fig. 2. — Three Gaussian fits to the optical depth profile in the direction of 3C123 demonstrating the fit improvement with additional 
components. Top p anels : absorption spectrum (exp(—r)) with Gaussian components overlaid; Bottom panels: fit residuals with ±ov(u) 
overlaid (see Section |3. 2 D . (a): base fit, 5 components, (b): additional component added (thick red dashed), 6 total, improving the fit with 
99% confidence by the F-test, so the component is retained, (c): additional component added (thick red dashed), 7 total, improving the fit 
with 61% confidence by the F-test, so the additional component is rejected. 


To compare the Xred values of different possible Gaus¬ 
sian fits, such as those shown for 3C123 in Figure (2j we 
employ an F-test (e.g., Westmoquette et al. 2007, Daw¬ 
son et al. 2011) using the IDL procedure mpftest.prcQ. 
The F-test allows us to quantify the improvement in 
Xred as a function of added DOF to the Gaussian model. 
First, the procedure computes the F-value, which is de¬ 
fined as a ratio of Xred an( l v va l ues f° r two models: 


Xred.lM 

XL,2/^2’ 


( 6 ) 


where the subscript “1” denotes the base model which in¬ 
volves the fewest number of Gaussian components, and 
“2” denotes a model including an additional component 
(i.e., three additional DOF). The probability density 
function pf(f\v 1 ,^ 2 ) for all F-values (/) described by 
DOF and v 2 is given b}0, 


1 ,^ 2 ) = 


r[(iq + V 2 /2)\ ( V 1 


i/l/2 


— X- 


jl/2(i/ 1 -2) 


F(^i/2)r(t' 2 /2) \v 2 ) " (1 + /Vi/V 2 ) 1/2( 

(7) 

From this we calculate the probability that a value / 
drawn from p/(f', v\, v 2 ) equals or exceeds the given value 
of F (Equation [6]), which is given by, 


pOO 

P F { F;i/ 1 ,v 2 ) = df Pf(f',vi,v 2 ), 

J F 


( 8 ) 


where Pf (or confidence level Cf = 1 — Pf) is a mea¬ 
sure of the probability that the added component merits 
inclusion in the model. 

For the majority of sources, the inclusion of addi¬ 
tional Gaussian components results in Cf > 0.99 be¬ 
fore C F drops far below Cf ~ 0.9 and indicates low 


5 http://purl.com/net/mpfit 

6 NIST/SEMATECH c-Handbook of Statistical Methods, 
http://www.itl.nist.gov/div898/handbook/ , 01/2015. 


confidence in the most recent component-addition. In 
three cases (sources 3C111, 3C133 and 3C459), the fit 
parameters varied strongly with the initial guesses, and 
we observed that the fits remained stable between the 
final component-additions if we applied a cutoff value of 
Cf = 0.97. Therefore, to apply a uniform standard to 
all sources, we use a strict threshold of Cf = 0.97 for 
including components, noting that in all but three cases, 
the effective cutoff is Cf > 0.99. 

As an example, in the case of 3C123, the confidence 
level of including one additional component (Figure [2b) 
to the base fit (Figure [2b) is C F = 0.99, indicating 
that the additional component exists with high confi¬ 
dence. The confidence level of including a second addi¬ 
tional component (Figure [2jc) in addition to the new fit 
(Figure Eb) is Cf = 0.61, and therefore we reject this 
additional component. 

3.4. Component Selection: Emission 

vi +^ he selection of Gaussian components for the expected 
emission profiles is more difficult. Taking the absorption 
line fits as fixed in central velocity, peak optical depth 
and line width (FWHM), we simultaneously solve for the 
Gaussian parameters of additional emission components 
and the spin temperatures of the absorption components 
given emission information along the line of sight (Equa¬ 
tions [3] and Considering the effects of stray radia¬ 
tion on the emission profiles and the complexity of the 
fits, we select a minimum number of components that 
brings the fit residuals below the 3aT exp (v) level. We 
then apply a similar F-test as in the absorption compo¬ 
nent selection and add additional components until the 
confidence level of the fit improvement drops below the 
cutoff (Cf = 0.97). 

Given the final selection of Gaussian components, we 
then solve the radiative transfer equations again for all 
possible orderings of the N absorption components along 
the line of sight. We repeat this line of sight iteration for 
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Fig. 3.— Gaussian fitting decomposition for 4/31 completed sources. The sources are displayed in two-panel plots, where the top panel 
displays the emission profile (from HT03, a2770, or the LAB survey) and the bottom panel displays the VLA absorption profile. The fits 
are overlaid in several forms, including total CNM and WNM contribution and individual components, according to the inset legend in 
panel (a). Residuals in the fits are offset at the bottom of each panel, with ±o’j' exp (?i) (top, see Section 2.4) and ±oy(u) (bottom, see 
Section 3.2) to illustrate the goodness of fit. 
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Fig. 4.— Zoom-in of a broad, warm component (green) in 
the complex profile towards PKS0531+19- Velocity: t>o,LSR = 
—24.96 it 4.25kms -1 , FWHM: An = 11.5 it 9.3kms -1 , peak op¬ 
tical depth: r pea k = 0.15 it 0.17, maximum kinetic temperature: 

max = 2900 it 5000 K, spin temperature: T s = 1451 it 263 K. 
The dotted line indicates a T (v). 

all permutations of Fk , the fraction of WNM unabsorbed 
by intervening CNM (see Section 3.1) among the I\ emis¬ 
sion components. If the total number of absorption com¬ 
ponents was large (defined here as greater than or equal 
to 6), we only varied the order of the most blended com¬ 
ponents at the center of the spectrum. The effect of 
different CNM cloud order on T s is very small compared 
with the effect of different values of Fk per emission com¬ 
ponent, therefore we held the more isolated components 
in these complex spectra fixed in their LOS order with 
little consequence to the derived T s values for the overall 
fit. 

4. SELECTED EXAMPLES 

To illustrate the fitting method further, we select four 
example sources and display their spectra and Gaus¬ 
sian decompositions in Figure [3] The top panel of each 
plot contains the emission profile, T exp (u), from HT03, 
a2770, or LAB (thin solid), the total fit to T exp (v) (thick 
solid), the total absorption contribution to the emission 
profile, Tb,abs(^) (Equation 1, dashed), and the total 
emission-only contribution, Tb,em(iO (Equation 1, triple 
dot-dash). The bottom panel of each plot contains the 
VLA absorption profile (thin solid), total fit to the ab¬ 
sorption profile (thick solid), and individual Gaussian 
components (dashed). The residuals of each fit are plot¬ 
ted at the bottom of each panel, with noise envelopes 
a T exp (v) (Section 2.4) and ±a T (v) (Section 3.2) to illus¬ 
trate the goodness of fit (dotted). 

We compare in detail our Gaussian decomposition with 
overlapping Hi absorption line surve ys, including HT03, 
Mohan et al. (2004) and iRov et all (l2013bO . All three 
surveys fit Gaussian functions to their absor ption l ines 
a nd f ollowed different strategies. For example. IRov et all 
(12013b!) fit increasing numbers of Gaussian functions un¬ 


til Xied f° r the fit was as close to 1 as possible. This 
resulted in a minimum of 1, maximum of 20 and of 
median of 7 Gaussian components over their 30 spec¬ 
tra. This approach requires detailed knowledge of the 
spectral noise properties, as well as the assumption that 
Gaussian functions are the proper representation of all 
components, which can fail if line broadening is not due 
to thermal and/or turbulence motions. Our method also 
makes these assumptions, but we fit fewer components 
per spectrum (minimum of 1, maximum of 11 and me¬ 
dian of 5) to not “over-fit” at the expense of inferior y/ ed 
values for the fits. In another approach, HT03 stress 
that human judgment is necessary for producing reliable 
Gaussian fits. For example, two fitted components with 
similar central velocities and FWHMs could in reality be 
a single component with a non-Gaussian intrinsic line- 
shape. No Gaussian decomposition is a unique solution, 
and it is impossible to prove which solution is most phys¬ 
ically relevant, but the differences between our respective 
Gaussian fitting strategies will affect the comparison of 
our results. We emphasize that our F-test method for 
component selection imposes a uniform standard to a 
non-unique fitting problem, and aims to make our deci¬ 
sions reproducible. 

4.1. 3C138 (Figure Hi) 

3C138 is included in the Millennium Arecibo survey 
(HT03), and they identified many similar, strong com¬ 
ponents, as we do, in their Gaussian decomposition. As 
seen in Figure [3Jr, in addition to narrow cold compo¬ 
nents at velocities 9.1, 6.4, 1.6 and -0.5 km/s, we iden¬ 
tify a shallow broad underlying component, centered at 
1.9 ± 0.3 km/s with a FWHM of 15.1 ± 0.6kms -1 . 
This corresponds to a maximum kinetic temperature of 
Tk ,max = 5000 K, and with a fit to the emission pro¬ 
file (from HT03), corresponds to a spin temperature of 
456 ± 40 K. The maximum spin temperature derived by 
HT03 for an absorption component in the direction of 
this source is T s = 380±23K, corresponding to a similar 
component at 1.8kms -1 . 

This source is also included in the Roy et al. (2013b) 
WSRT/GMRT/ATCA survey. Whereas we fit a total 
of 6 components before the fit ceases to improve signif¬ 
icantly by the F-test, Roy et al. (2013b) fit 13 compo¬ 
nents. Given the complex nature of the spectrum, and 
the large difference in total numbers of components, there 
are few closely matching components between the two 
decompositions. 

4.2. 3C286 (Figured) 

With a sensitivity of cr T = 9 x 10 -3 , HT03 were un¬ 
able to detect Hi absorption in the direction of 3C286. 
Stanimirovic & Heiles (2005) demonstrate that with ad¬ 
ditional integration times at the Arecibo Observatory, 
weak absorption lines below the HT03 sensitivity thresh¬ 
old are easily detected, which can have important im¬ 
plications for HT03’s measured fractions of Hi in the 
CNM, WNM and thermally unstable phases. For 3C286, 
they find 3 velocity components at —28.8, —14.3 and 
—7.4kms -1 , with ax = 106, 115 and 315K and 
T s = 89 ± 7, 37 ± 4 and 30 ± 20 K. We find similar veloc¬ 
ity components, at —28.5, —14.3 and —7.3kms -1 , with 
Xfc jlnax = 116, 220 and 405 K and T s = 63±13, 84±12 and 
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27 ± 9K. Roy et al. (2013b) find the same three compo¬ 
nents, in addition to a broad line centered at — 8.4 km s" 1 
with peak optical depth r = 0.001 and Tfc jmax = 5140 K. 
Their achieved sensitivity of ay = 3 x 10 -4 allows them 
to detect this weak component whereas our sensitivity of 
<j T = 7 x 10 -4 puts it below our detection limit. We note 
that they do not discuss the effect of bandpass instability 
on their detection of broad, weak lines. 

4.3. 4C16.09 (Figure^) 

Although 4C16.09 is not included in the Millennium 
survey (HT03), it is included in the Mohan et al. (2004) 
GMRT survey of Hi absorption. They achieve a mean 
survey sensitivity in optical depth ofcr T ~3xlO -3 and 
the spectra have velocity resolution of 2.1 kms -1 over 
roughly 350 km s -1 . Our absorption profile agrees well in 
shape with their results. They detect two narrow absorp¬ 
tion components, at velocities 7.5 and 0.2 kms -1 with 
corresponding spin temperatures T s = 30 K and 45 K re¬ 
spectively (Mohan et al. 2004). We detect components 
at similar velocities, 7.6 and 0.7kms , in addition to 
four additional absorption components (see Figure [3J;). 
We identify two additional, broad emission components 
which are prominent in the expected emission profile and 
produce a maximum spin temperature of T s = 240±10 K. 
In comparison, Roy et al. (2013b) have the same sensi¬ 
tivity of ay = 8 x 10 -4 and find 8 Gaussian components. 

4.4. PKS0531+19 (FigureW) 

The profile in the direction of PKS0531+19 is highly 
complex. Detected by HT03 in the Millennium Survey, 
they identify four narrow absorption components at ve¬ 
locities 9.6, 5.5, 1.8 and —23.1 kms -1 (HT03). We detect 
features at similar velocities, in addition to two broad, 
shallow components centered at 21 and — 25 kms -1 cor¬ 
responding to spin temperatures of T s = 1129 ± 74 K and 
1451 ± 263 K. In Figure Q] we zoom in on the second of 
these components, shown by the thick green line. This 
source is a good example of the benefit of high sensitivity. 
The RMS noise in the absorption profile is ay = 5 x 10 -4 , 
so that the broad component with peak optical depth 
equal to r = 0.0015 is a 3er detection. HT03 achieve a 
sensitivity of ay ~ 0.0015 in their absorption profile to¬ 
wards this source, and therefore did not fit a Gaussian 
function to this line given that it fell below their ler level. 

5. PROPERTIES OF COLD AND WARM NEUTRAL GAS 

Following complete decomposition of all absorption 
and emission profiles with Gaussian functions, we can 
use estimates of temperature and column density to in¬ 
fer properties of the CNM and WNM. 

5.1. Temperature 

The FWHM (Av) of each Gaussian function deter¬ 
mines its maximum kinetic temperature, T^ max , by 

T k , max = m H /{8k B ln(2)) x Av 2 = 21.866 x Av 2 , (9) 

for hydrogen mass m B , Boltzmann’s constant k B , and 
FWHM Av in kms -1 (e.g., HT03). This quantity con¬ 
tains contributions from both thermal broadening and 
turbulent gas motions so that, 

r fc ,max=Tfc+ (10) 

2 k B 


where Tk is the true kinetic temperature of the gas (e.g. 
Liszt 2001). Therefore, X/^max is an upper limit to the 
kinetic temperature in the presence of turbulent broaden¬ 
ing. As shown by Heiles & Troland (2003b), nonthermal 
motions in the CNM are characterized by a turbulent 
Mach number (M t ) of M t ~ 3, implying that CNM gas 
is supersonic. To estimate M t , we solve Equation m 
for u 2 urb , and then assume that this line of sight tur¬ 
bulent velocity is \/(3) times the three-dimensional tur¬ 
bulent velocity (vturb, 3 D)- The turbulent Mach number 
M t is given by ^turb 3 D/^s> where C s is the isothermal 
sound speed. We use C s because thermal equilibrium 
is established quickly in the CNM, which dominates the 
absorption profiles (Heiles & Troland 2003b). Following 
Heiles & Troland (2003b), we adopt C 2 = kT s /lAm B 
for mean atomic weight 1.4mn. Therefore, the turbulent 
Mach number is given by, 

= = 4.2(%^-lV (ii) 

^ s V x s / 

which is Equation (17) from Heiles & Troland (2003b). 
In Figure [Sj we display a histogram of M t for all 
21-SPONGE absorption-detected Gaussian components, 
which has a peak at M t ~ 3 (median M t = 2.86) and 
agrees with the results from HT03. The distribution is 
strongly peaked, and 40% of all components have Mach 
numbers between 2 and 4. 

The spin temperature is determined by the ambient 
radiation field, collisions between hydrogen atoms, elec¬ 
trons and protons, and Lya scattering. For the CNM, 
high densities allow the collisional contribution to domi¬ 
nate 21-cm excitation so that T s is expected to be equal 
to Tfc. At low densities, i.e. in more diffuse WNM gas, 
collisions are not sufficient to thermalize the 21-cm tran¬ 
sition and therefore T s is expected to be less than Tk (e.g. 
Liszt 2001). 

In the case of 12 sources, one or more of the Gaussian 
functions fit to the absorption spectrum did not have a 
corresponding spectral feature at the same velocity in 
emission (likely due to the cold component occupying a 
very small solid angle, Heiles & Troland 2003a). There¬ 
fore, the converged T s values for these components are 
unreasonably low, < 1K. Given the uncertain presence 
of stray radiation in the Arecibo emission profiles, we 
take a conservative approach for estimating T s for these 
components. We adopt the median M t value measured 
over all fitted Gaussian functions, M t = 2.86 and solve 
Equation (fill) for T s and its error given Tfc max and its 
error. 

5.2. Column Density 

We estimate the Hi column density of the absorbing 
Hi using the expression: 

N(Hi, ABS) =C 0 T S j t(v) dv (12) 

where Co = 1.823 x 10 18 cm -2 K -1 (kms -1 ) -1 , T s is the 
derived spin temperature and t(v) is the optical depth as 
a function of velocity in km s -1 . For each n th absorption 
component, we use parameters from the Gaussian fit to 
calculate the column density, so that, 

N(Hi, ABS) n = 1.064 xC 0 x T SiU x To, pea k,n x Av n , (13) 
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Mach number, M, 

Fig. 5. — Histogram of turbulent Mach number Mt for all 
abso rptio n-detected Gaussian components, calculated from Equa¬ 
tion <| 11 f) (Heiles & Troland 2003b Equation 17). 

where (T Si „,ro lP eak,n, Aii„) are the (spin temperature, 
peak optical depth and FWHM in km s _1 ) of the n th 
component and 1.064 converts the product to the area 
under a Gaussian function with the given height and 
width. For the emission components not detected in ab¬ 
sorption, we estimate the column density as, 

N(Hi, EM) = C 0 j T b {v) dv, (14) 

where T B is the brightness temperature. Using the Gaus¬ 
sian fit to the emission profile, we estimate the column 
density associated with each k th non-absorbing compo¬ 
nent as, 

N(Hi, EM) k = 1.064 x Co x T 0 ,fc x Av k , (15) 

where To.fc and Av k are the peak and FWHM in km s -1 
of the k th component. 

6. COMPARISON WITH PREVIOUS SURVEYS 

In Figure [HJ for all absorption-detected Gaussian func¬ 
tions we display histograms of: (a) spin temperature 
(T s ), (b) maximum kinetic temperature (T^ max ), (c) 
peak optical depth (To, pea k), (d) column density for 
all CNM clouds, defined as those with T s < 200 K, 
(N(Hi,CNM)), and finally (e) column density for all 
warm clouds, defined as those emission-only components 
and absorption-detected components with T s > 200 K 
(N(Hi, warm)). We overlay the same quantities for the 
79 sources from the Millennium survey (HT03) in the 
panels of Figure [6] for comparison. 

The superior sensitivity of 21-SPONGE allows us to 
detect absorption lines with smaller r pea k and/or higher 
T s than in HT03, which explains the poor agreement be¬ 
tween the 21-SPONGE and HT03 distributions in Fig¬ 
ures Et and El:. We are therefore more sensitive to low- 
N(Hi, CNM), as shown in FigureEH and in Equation (full) . 

In Figure Et, we detect a strong tail of high-T s com¬ 
ponents in comparison with HT03. This results in our 
median T s = 77 K, which is higher than their median 
value of T s = 48 K. The maximum T s we detect is 


T s = 1451 ± 263 K, compared with the maximum value 
measured by HT03 of T a = 656 K, demonstrating that 
it is possible to directly measure T s from warm gas with 
improved sensitivity in Hi absorption. 

However, although we have the sensitivity to detect 
Hi at higher spin temperatures than 1500 K directly in 
absorption, we only detect 32/157 Gaussian components 
above the expected range for the CNM ( T s = 40 — 200 K 
from Wolfire et al. 2003), and only 4/157 above T s = 
1000 K. This is a surprising result from 21-SPONGE so 
far, given that HT03 estimate (based on indirect mea¬ 
surements from Hi emission only) that 48% of com¬ 
ponents should be in the thermally unstable range of 
500 < T s < 5000 K. Our current estimate for the ther¬ 
mally unstable fraction in this same range- ~ 20% by 
number- is more consistent with estimates from collid¬ 
ing flow models (e.g., Audit & Hennebelle 2005). 

In Figure [Tj we present cumulative distribution func¬ 
tions of the data in Figure Et- A K-S test indicates that 
the two distributions are highly unlikely to be drawn 
from the same parent population (K-S = 0.23, p = 0.00). 
The 21-SPONGE distibrution is much shallower at high 
temperature, indicating again that we are detecting more 
gas at higher T s than was done by HT03 using similar 
analysis methods. 

In Figures El 3 and Et we also include Gaussian de¬ 
composition results from the WSRT/GMRT/ATCA sur¬ 
vey presented in Roy et al. (2013b) (grey shaded his¬ 
tograms; see Table 1 for survey information). Our two 
surveys have similar sensitivity in absorption. Follow¬ 
ing Gaussian decomposition of their absorption spectra, 
Roy et al. (2013b) do not estimate T s , AT(Hi, CNM) or 
N(Hi, warm) for their Gaussian components. As dis¬ 
cussed in Section [U Roy et al. (2013b) fit many more 
components per spectrum than we do, and therefore they 
detect many more low-r lines than either we or HT03 do, 
as seen by the strong low-r tail in the gray-shaded his¬ 
togram in Figure Et- 

We agree well with both Roy et al. (2013b) and HT03 
in T/^max, shown in Figure Et>- From Equation (El), this 
implies that we are sensitive to similar absorption line 
widths. We find a median value of Tfc max = 160 K, com¬ 
pared with HT03’s median value of T ktI nax = 220 K and 
Roy et al. (2013b)’s value of Tfc imax = 160 K. We find 
that 5.7% of components have Tfe.max < 40 K, which is 
below the minimum theoretically expected kinetic tem¬ 
perature for the CNM (e.g. Wolfire et al. 2003), com¬ 
pared wtlr 4.7% by HT03 and 7.0% by Roy et al. (2013). 
At the extreme end of the low-Tfc jlnax values, we find one 
weak component with T kjnax = 9 ± 4K in the direc¬ 
tion of 3C345. These low-temperature components can 
be explained by ineffective photoelectric heating from 
dust and polycyclic aromatic hydrocarbons (PAHs) along 
those sightlines, which allows for the existence of very 
cold clouds with T kjiaaK ~ 10 K (Wolfire et al. 1995, Heiles 
& Troland 2003b, Roy et al. 2013b). 

In Figures EH and Et, it is clear that we generally 
agree with HT03 in measurements of N(Hi, CNM) and 
N(Hl,warm). However, our superior sensitivity in op¬ 
tical depth allows us to detect lines with much lower 
N(Hi, CNM). For example, we are able to detect an in¬ 
dividual absorption line with N(Hi, CNM) = 3 ± 1 x 
10 16 cm~ 2 , whereas HT03 measure a minimum column 
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Fig. 6. — Histograms of various properties from the Gaussian fits to the current 31 sources in 21-SPONGE (thick black), including: (a) 
T s (K), (b) Tk imax (K), (c) peak optical depth in absorption, ro )Pe ak 5 (d) N(Hi, CNM): column density of all CNM components, defined as 
those with T s < 200 K, (d) N(Hi,warm): column density of all “warm” components, defined as those with T s > 200 K and those detected 
only in emission. The results for the Millennium Survey (HT03, thin black) and the WSRT/GMRT/ATCA survey by Roy et al. (2013b) 
(shaded gray) are included for comparison. Results from Roy et al. (2013b) are not included in panels a, d and e because they do not 
compute these parameters for individual Gaussian components. 


density of 2 x 10 18 cm -2 for individual CNM clouds. An 
important result from the final 21-SPONGE data re¬ 
lease will be constraints on the low-N(Hi) end of the 
CNM column density distribution. McKee & Ostriker 
(1977) predicted a lower-limit to CNM column density 
of N(Hi, CNM) ~ 10 17 — 10 18 cm -2 , derived from bal¬ 
ancing radiative losses and conductive heating at the in¬ 
terface between the CNM and WNM at various possi¬ 
ble temperatures. As discussed by Stanimirovic & Heiles 
(2005), observational constraints on this lower limit (and 
its prevalence) can distinguish between different forma¬ 
tion and evolution histories of these structures, whether 
by evaporation into the surrounding medium, or tran¬ 


sient, dynamic formation at shock interfaces, or conden¬ 
sation from turbulent flows. 

6.1. Line of sight properties 

In Table [3] we compute line of sight column density- 
related quantities for the current 21-SPONGE sources, 
where the subscript “20” indicates units of 10 2O cm” 2 . 
These quantities include: (a) the total CNM column 
density (EN(Hi, CNM) 2 o), equal to the sum of the col¬ 
umn densities of the absorption-detected Gaussian com¬ 
ponents with T s < 200 K; (b) the total WNM col¬ 
umn density (EN(Hi, warm) 2 o), equal to the sum of the 
column densities of the absorption-detected Gaussian 
components with T s > 200 K plus the column densi- 
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TABLE 3 

Line of sight column density parameters 


Source 

EN(Hi,CNM) 2 o 

EN(Hi, warm) 20 

N(Hi, total) 2 o 

/cnm 


(10 20 cm -2 ) 

(10 20 cm“ 2 ) 

(10 20 cm -2 ) 



4C32.44 

0.17 

± 

0.09 

1.19 

± 

0.03 

1.36 

=t 

0.09 

0.12 

=t 

0.06 

3C286 

0.06 

± 

0.02 

1.28 

± 

0.09 

1.34 

± 

0.09 

0.05 

=t 

0.02 

4C12.50 

0.36 

± 

0.11 

2.01 

± 

0.03 

2.37 

± 

0.11 

0.15 

=t 

0.05 

3C273 

0.05 

± 

0.19 

2.32 

± 

0.03 

2.37 

± 

0.19 

0.02 

± 

0.08 

3C298 

0.06 

± 

0.37 

2.33 

± 

0.09 

2.39 

± 

0.38 

0.03 

± 

0.15 

4C04.51 

0.24 

± 

0.06 

4.31 

± 

0.11 

4.55 

± 

0.13 

0.05 

± 

0.01 

3C237 

0.17 

± 

0.20 

2.43 

± 

0.11 

2.60 

± 

0.23 

0.07 

± 

0.08 

3C225A 

1.19 

± 

0.30 

3.05 

± 

0.06 

4.24 

± 

0.30 

0.28 

± 

0.07 

3C225B 

0.98 

± 

0.18 

3.57 

± 

0.11 

4.55 

± 

0.21 

0.22 

± 

0.04 

3C345 a 

0.00 

± 

0.00 

1.01 

± 

0.13 

1.01 

± 

0.13 

0.00 

± 

0.00 

3C327.1 

1.76 

± 

1.17 

6.96 

± 

0.11 

8.72 

± 

1.18 

0.20 

± 

0.14 

3C147 

7.44 

± 

1.79 

13.15 

± 

0.78 

20.59 

± 

1.95 

0.36 

± 

0.09 

4C33.48 

1.54 

± 

3.00 

14.22 

± 

0.48 

15.76 

± 

3.04 

0.10 

± 

0.19 

3C154 

9.75 ± 14.94 

36.19 

± 

0.49 

45.94 ± 14.94 

0.21 

± 

0.33 

3C410 

30.39 

± 

6.64 

29.29 

± 

0.81 

59.69 

± 

6.69 

0.51 

± 

0.13 

B2050 

3.00 

± 

2.02 

30.86 

± 

1.33 

33.86 

± 

2.42 

0.09 

± 

0.06 

3C409 

5.33 

± 

3.49 

26.36 

± 

0.55 

31.69 

± 

3.53 

0.17 

± 

0.11 

PKS0531 

8.39 

± 

9.61 

20.45 

± 

1.48 

28.84 

± 

9.72 

0.29 

± 

0.35 

3C111 

3.99 

± 

3.78 

32.44 

± 

2.29 

36.43 

± 

4.42 

0.11 

± 

0.10 

3C133 

9.76 

± 

4.02 

26.22 

± 

0.46 

35.98 

± 

4.05 

0.27 

± 

0.12 

3C138 

6.01 

± 

2.83 

20.42 

± 

0.58 

26.43 

± 

2.89 

0.23 

± 

0.11 

3C123 

5.15 

± 

4.41 

19.41 

± 

1.90 

24.55 

± 

4.80 

0.21 

± 

0.18 

3C433 

1.15 

± 

2.09 

9.23 

± 

0.13 

10.38 

± 

2.10 

0.11 

± 

0.20 

3C120 

10.86 

± 

2.24 

11.36 

± 

0.26 

22.22 

± 

2.26 

0.49 

± 

0.11 

3C48 

0.70 

± 

0.36 

4.25 

± 

0.19 

4.95 

± 

0.41 

0.14 

± 

0.07 

4C16.09 

3.27 

± 

1.10 

8.83 

± 

0.24 

12.10 

± 

1.13 

0.27 

± 

0.09 

3C454.3 

1.99 ± 15.74 

5.92 

± 

0.63 

7.90 ± 15.75 

0.25 

± 

2.05 

J2232 

1.21 

± 

0.37 

4.29 

± 

0.17 

5.50 

± 

0.41 

0.22 

± 

0.07 

3C78 

3.81 

± 

1.91 

9.93 

± 

0.27 

13.74 

± 

1.93 

0.28 

± 

0.14 

3C459 

1.17 

± 

0.49 

5.12 

± 

0.04 

6.29 

± 

0.49 

0.19 

± 

0.08 


IV(Hi, CNM) 20 = 3 ± 1 x 10 16 cm" 2 



Fig. 7.— Cumulative distribution function of the T a values 
shown in Figure [6^ for 21-SPONGE and HT03. 

ties of all emission-only Gaussian components (which 
is a combination of both thermally-stable WNM and 
thermally-unstable gas); (c) the total Hi column den¬ 
sity per iine of sight (N(Hi, total) 2 o = £N(Hi, CNM) 2 o + 
EN(Hi,warm) 2 o); (d) the CNM fraction, /cnm, defined 


Fig. 8.— “Box and whisker” plot comparing the CNM fractions 
(/cnm) for 21-SPONGE with HT03 and Roy et al. (2013). Each 
box spans the 25 th through 75 th percentiles of each distirbution, 
the median value is indicated by the horizontal line within each 
box, and the bars extend to the minimum and maximum values. 


as: 
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_ £N(Hi, CNM) 
/CNM ~ N(Hi, total) ' 


(16) 


To compare our results, we obtained or calculated the 
same quanitites for the HT03 and Roy et al. (2013b) 
samples. For HT03, we used their published fit in¬ 
formation to calculate all quantities in the same man¬ 
ner as for 21-SPON GE, described abo ve. To calcu¬ 
late £N(Hi, CNM) 20 , IRov et all (l2013blf began by as¬ 
suming that all absorbing gas is from the CNM with 
T s = 200 K and solved Equation m to produce upper 
limits to £N(Hi, CNM) 2 o for each line of sight. We es¬ 
timated £N(Hl, warm) 2 o from their results by subtract¬ 
ing £N(Hl, CNM) 2 o from 1 V(Hi, total) 2 o, which they cal¬ 
culated using the “isothermal” method (Spitzer 1978, 
Dickey & Benson 1982, Chengalur et al. 2013). 

We gene r ally ag ree well with the results from HT03 and 
IRov et ahl (j2013bi l. All three studies measure the same 
range in £N(Hl, CNM) 20 , between ~ 0.01 and ~ 10. We 
find a median £N(Hi, CNM) 2 o = 1.76, compared with 
1.45 by HT03 and 1.76 by Roy et al. (2013). From Fig¬ 
ure [U we detect many more components in absorption 
with T s > 200 K, which have high 1 V(Hi, ABS) n accord¬ 
ing to Equation m and which will increase the total 
£N(Hi, warm) 2 o along the line of sight. We measure me¬ 
dian £N(Hi, warm) 2 o = 8.83, which is much higher than 
HT03’s median value of 3.97 and Roy et al. (2013)’s me¬ 
dian value of 1.41. We measure a median total column 
density per LOS of lV(Hl, total) 2 o = 10.4 compared with 
5.54 by HT03 and 4.33 by Roy et al. (2013). 

Figure[5]displays a “box and whisker” plots of /cnm for 
21-SPONGE, HT03 and Roy et al. (2013a). By inspec¬ 
tion, our results are consistent with the /cnm distribu¬ 
tion found by HT03 and Roy et al. (2013a). HT03 do not 
detect Hi absorption in the direction of 15/79 sources, 
resulting in many values of /cnm = 0. Whereas we 
only have one non-detection, we find 9/31 sources with 
/cnm < 0.1. We measure median /cnm = 0.20 com- 
paredwith 0.23 by HT03. We appear to find lower /cnm 
than IRov et al.1 (l201 3bl. Out of our 31 sources, we find a 
maximum /cnm ~ 0.51, whereas IRov et all (l2013bll find 
much higher values, all the way up to /cnm ~ 1 with 
median 0.51. As discussed above, their £N(Hi,CNM) 20 
are strong upper limits, relying on the assumption that 
all absorbing gas has T s = 200 K, which neglects gas 
along the line of sight at higher and lower temperatures. 
Therefore, their /cnm values, calculated by dividing the 
upper limits to £N(Hi,CNM) 2 o by the total isothermal 
column density (Chengalur et al. 2013), are all strong up¬ 
per limits. 


7. COMPARISON WITH SYNTHETIC OBSERVATIONS 

Recently, Kim et al. (2014) used a hydrodynamical sim¬ 
ulation of a Milky Way-like Galactic disk to construct 
synthetic Hi absorption and emission profiles for thou¬ 
sands of sight lines probing the turbulent, multiphase 
ISM. They calculate “observed” T s , 7/, N(Hi) and /cnm 
for all lines of sight and found that the majority of sight 
lines with average LOS spin temperature below 200 K 
have a CNM fraction between 0.4 and 0.7 (98% have 
/cnm < 0.7; C-.G. Kim et al. 2014 private communica¬ 
tion). Whereas we measure several CNM fractions within 
this range, by inspection of Figure [8] and Tabled we have 


more measurements below /cnm = 0.3 than predicted 
by the Kim et al. (2014) results. However, in agreement 
with HT03, Kim et al. (2014) also found many WNM- 
dominated sightlines with /cnm = 0. 

Interestingly, Kim et al. (2014) found that 
synthetically-observed line-of-sight (harmonic mean) 
T s values agree within a factor of 1.5 with “true” 
harmonic mean T s values computed from the simu¬ 
lation over the full range of integrated optical depth 
(10 -3 < frdv < 10 2 kms _1 ). They argue that this 
implies that there is little overlap between CNM clouds 
along a line of sight within their lkm s _1 velocity 
channels. However, by inspection of our absorption 
profiles and fits (e.g. Figure 4), it is clear that cloud 
complexes can overlap significantly along every real 
observed line of sight. 

In Figure [9] we plot CNM fraction versus the opti¬ 
cal depth-weighted LOS spin temperature, T Si0 b s - This 
weighted spin temperature, T SjC1 b s , is calculated from 
Equation (15) of Kim et al. (2014), which was derived 
from the “isothermal” column density estimate (Spitzer 
1978, Dickey & Benson 1982, Chengalur et al. 2013, Kim 
et al. 2014), and is given by, 


Ts,obs — 


f T ( V ) i.^(> 

J r{v)dv 


(17) 


We display the error in the weighted T s values when 
they exceed the symbol size. The squares represent 
same quanitites with a cutoff for the CNM fraction at 
T s = 350 K instead of T s = 200 K. This increases the 
value of /cnm in only a few cases, because we only detect 
6% of all components in the 200 <T S < 350 K regime by 
number. This shows that our choice of temperature cut¬ 
off does not strongly affect the observed trends in /cnm- 

Kim et al. (2014) plotted the same quantities (their 
Figures 7a and 8a) and found that /cnm, defined by a 
7/ cutoff at 184 K, varies with T Si0 b s roughly as /cnm ~ 
T C /T S)0 b s for 50 < T c < 100 K. These T C /T S)0 b s curves 
are shown in Figure [9] in dashed purple for T c = 20, 50 
and 100 K. Although we see a similar tight trend be¬ 
tween CNM fraction and T Sj0 b s , we find points below 
the 50/T SjO bs curve, indicating again that we are measur¬ 
ing lower /cnm and motivating us to extend the dashed 
curves down to 20 /T s . Kim et al. (2014) possibly pro¬ 
duced higher /cnm because they did not include Hi-H 2 
chemistry, and therefore could have overestimated /cnm 
if some cold, dense Hi should have transitioned to H 2 
form. This was also noticed by Stanimirovic et al. (2014) 
in an Hi emission and absorption line study of the Perseus 
molecular cloud at the Arecibo Observatory. 

Overall, we generally agree with the trend between 
CNM fraction and optical depth-weighted LOS T s 
{T s , obs) produced by the synthetic observations of the 
Kim et al. (2014) simulation, although our observations 
increase the scatter in the trend. The agreement is 
very encouraging considering that our CNM fractions 
rely on temperature calculation by Gaussian decompo¬ 
sition, whereas theirs are taken directly from their simu¬ 
lated lines of sight. Kim et al. (2014) derived an expres¬ 
sion describing this trend based on several assumptions: 
T a = Tfc in the CNM, the “isothermal” column density 
equals the true column density (truest for low integrated 
r sightlines), and the synthetically observed harmonic 
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mean T a equals the true harmonic mean T s (where the 
agreement is best for r < 0.1, above which there is scatter 
of up to 1.5). Finally, in computing /cnm, they assume 
that the CNM fraction is not “extremely small” in order 
to express the /cnm as a ratio of temperatures instead 
of column densities (see their Section 3.2, Equation (16), 
Kim et al. 2014). This last assumption may be the most 
important to increasing the trend scatter, because we 
(and HT03) find many sightlines with both small /cnm 
and small T Si0 b s - 

The differences between our results are likely due in 
part to the fact that Kim et al. (2014) do not include 
magnetic fields or chemistry to map the transitions be¬ 
tween Hi and H 2 , but the similarities found so far suggest 
that only small adjustments to the simulations may be 
necessary to match the observed /cnm- We also empha¬ 
size that the spatial resolution of a simulation is key for 
resolving the CNM and appropriately calculating quan¬ 
tities like the CNM fraction (e.g. Audit & Hennebellc 
2005). In the future we will apply our Gaussian fitting 
methods to large volumes of synthetic data from Kim et 
al. (2014) and other higher spatial-resolution numerical 
simulations using Autonomous Gaussian Decomposition 
(Lindner et al. 2015). This will be essential to assessing 
the effect of our analysis techniques on these results and 
also understanding the biases of the numerical models. 

8. COMPARISON OF METHODS FOR ESTIMATING SPIN 
TEMPERATURE FROM HI ABSORPTION AND 
EMISSION SPECTRA 

Several independent methods are commonly used in 
the literature to measure T s from Hi absorption and 
emission observations. The first is the method we use 
in this work, whereby we estimate a single T s value for 
each Gaussian component fitted to both absorption and 
emission spectra (i.e., solving Equation 1). This method 
relies on initial guesses for number of components and 
their shapes, but it is the only method based on a phys¬ 
ical model. As is generally true for the 21 cm line, a 
Gaussian model is a good representation of the spectral 
line because the damping wings are insignificant. Al¬ 
though this assumption fails in some cases, shown in this 
study by the components with 7/ jmax » T s (6% have 
7/,max > (T s + 3<t)), the method is convenient and pro¬ 
duces reasonable fits to both simple and complex profiles. 


8.1. One phase per channel T s 

Assuming that the gas in each velocity channel is at a 
single temperature (the “one-phase” assumption), divid¬ 
ing the brightness temperature spectrum (Tb(v)) by the 
absorption profile (1 — e~ T ^) produces a spin tempera¬ 
ture spectrum, T s ^o{v ), via the following equation, 


T s ,o{v) 


T b (v) 

1 - e~ T 0) ‘ 


(18) 


From this it is possible to measure a single T s value for 
every velocity channel. For example, Roy et al. (2013a) 
use this method to compute T Sj o for each velocity chan¬ 
nel above 3<r significance in absorption at 1 kms -1 reso¬ 
lution. 


8.2. Line of sight harmonic mean T s 


Assuming that there is one value of T s along a full line 
of sight, the “LOS” harmonic-mean temperature, T/l 
from the absorption and emission profiles is given by, 


J T B [v)dv 
f (1 — e~ T 0))dv 


(19) 


For example, Kanekar et al. (2011) computed T s b for the 
same data from Roy et al. (2013) to distinguish between 
CNM and WNM-dominated lines of sight. 


8.3. The Slope Method 

It is also possible to measure T s for individual gas 
components by estimating the slope or “ridge line” of 
each linear feature in Tb{v) vs. 1 — e~ T ^ space (e.g. 
Mcbold et al. 1997, Dickey et al. 2000). The advantage 
of this approach compared with Gaussian decomposition 
is that the solution can be determined directly from the 
data without fitting iterations or initial guesses. How¬ 
ever, this method also requires the assumption that all 
gas in a given velocity range exists at a single tempera¬ 
ture, as is also assumed in the one phase per channel and 
LOS methods. To improve the slope method parameter 
estimation, Dickey et al. (2003) developed a two-phase 
nonlinear least squares fitting technique to determine the 
T s slopes, which we will call the “NLLSQ slope method”. 

Dickey et al. (2003) carried out a detailed comparison 
of T s estimation methods, including Gaussian decompo¬ 
sition and the NLLSQ slope method, and found that in 
cases of isolated, unblended CNM components, these two 
methods agree reasonably. In another comparison study, 
HT03 computed T s for a sample of their sources using 
a by-eye version of the slope method and similarly con¬ 
cluded that the two methods agree in the case of un¬ 
blended CNM components, even if there are blended, 
weaker WNM components present. Dickey et al. (2003) 
re-processed all 79 HT03 emission and absorption pairs 
using the NLLSQ slope method and found a median spin 
temperature per components of T s = 31 K, which is lower 
than the median value per fitted Gaussian functions by 
HT03 of T s = 48 K. This suggests that there may be 
systematic differences between the two methods. 

In the comparisons made by Dickey et al. (2003) and 
HT03, the T s values from the NLSSQ slope method tend 
to be lower than the Gaussian-derived T s . This is due to 
the fact that spectral features in Hi emission and absorp¬ 
tion profiles do not always have the same shape. If chan¬ 
nels in emission corresponding to narrow absorption have 
similar brightness temperature as channels with no cor¬ 
responding absorption, the slope of the opacity feature 
will be small. HT03 note that this is because component 
emission is determined by the product of both optical 
depth and T s , so that broad, warm components dom¬ 
inate emission and narrow, cold components dominate 
only in absorption. We will conduct a detailed compar¬ 
ison of the Gaussian and NLLSQ T s derivation methods 
upon completion of the 21-SPONGE data collection. In 
addtion, we plan to use synthetic spectra from numeri¬ 
cal simulations, where the input T s values are known, to 
assess and quanitfy their associated biases. 

8.4. Applying methods to 21-SPONGE data 

From a preliminary analysis of 5 21-SPONGE sources 
(3C286, 3C225A, 3C225B,“ 3C237, and 3C298), we find 
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Optical depth-weighted average T s per LOS 


Fig. 9.— CNM fraction (/cnm) versus optical depth weighted average T s per LOS (Equation I17[) . Purple dashed lines indicate the 
observed trends in CNM frac tion with LOS spin temperature from Kim et al. (2014) (their Figures 7a and 8a). Crosses indicate /cnm 
calculated from Equation Gii with a cutoff of T s = 200 K to define the CNM, squares indicate /cnm using a cutoff of T s = 350 K. Errors 
in /cnm are shown (solid vertical lines) when they exceed the symbol size. Dotted lines connect crosses and squares when the difference 
between the two values exceed the symbol size. 


that the NLLSQ slope method temperatures are mostly 
in the range 25 < T s < 100 K, which agrees well with 
the peak of our T s distribution (Figure [6ji). In Figure ITOl 
we compare the Gaussian fit, one-phase per channel, and 
LOS harmonic mean methods by applying them all to 
the 31 21-SPONGE emission/absorption spectral pairs. 
For the one-phase per channel method, for each source 
we used Equation (1151) and included the T s values for all 
channels with absorption signals above the 3er level. In 
Figure flOh. we plot histograms of T s for (1) the Gaus¬ 
sian decomposition method (thick green, “G”) (2) the 
one-phase per channel method (thin black, “O”) (3) the 
LOS harmonic mean method (shaded gray, “L”). In Fig¬ 
ure [TO}:> we display cumulative distribution functions for 
the three T s methods. 

Based on the CDFs of the three T s distributions in Fig¬ 
ure m>, the one-phase per channel and LOS-integrated 
methods agree well (K-S = 0.18,p = 0.35). Both meth¬ 
ods rely on the assumption that there is a single-phase 
medium in each measurement interval (either velocity 
channel or full line of sight), which over-simplifies the re¬ 
ality that the ISM is a complex, multi-phase mixture of 
gas at many temperatures. The one-phase assumption 
overestimates T s , because within each velocity channel 
(and line of sight) narrow components at lower tempera¬ 
tures can cumulatively produce a stronger absorption or 
emission signal, thereby biasing the average T s measure¬ 
ment to high values. As previously discussed by Dickey 
et al. (2000) and Dickey et al. (2003), the NLLSQ slope 


and Gaussian methods provide systematically lower val¬ 
ues of T s than the single-phase approximations. This is 
supported by strong disagreements in Figure H0lip=0. 00 
between both G and O, and G and L). 

The median values of the Gaussian, one-phase per 
channel and LOS methods are T S) q = 77 K, T St o = 340 K 
and T s ,l = 290 K. Wolfire et al. (2003) model the multi¬ 
phase ISM with detailed analytical prescriptions for heat¬ 
ing and cooling, and show that the CNM should have a 
spin temperature of 40 < T s < 200 K. The Gaussian fit 
method best reproduces this predicted range, although 
the T s values in Figure [TO] also contain components as¬ 
sociated with thermally unstable gas in addition to the 
dominant CNM. Furthermore, Figure [1(11 shows that the 
Gaussian fit method is the most successful in reproduc¬ 
ing the analytical prediction for the lower limit to CNM 
T s . The LOS method in particular misses a large fraction 
of the CNM T s distribution, as it predicts all gas should 
have T s > 200 K. As discussed by Dickey et al. (2009), 
the T s single-phase (i.e. O and L) T s methods provide es¬ 
timates of the mean temperature or the mixture of CNM 
and WNM along the line of sight, rather than physical 
temperatures of individual structures. 

9. SUMMARY 

In this paper we present methods, data and prelimi¬ 
nary results from the first half of 21-SPONGE, a large Hi 
absorption survey at the Karl G. Jansky VLA. The ab¬ 
sorption data are complemented by emission data from 
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T s (K) T s (K) 


Fig. 10. — Comparison of temperature derivation methods for the 21-SPONGE survey, (a): Histograms of T s fo r: G aussian decomposition 
of Hi absorption and emissio n li nes (thick green), “one-phase” T s per 0.42 km/s channel (thin black, Equation 1 1811 and integrated T s per 
LOS (shaded gray, Equation 1191 1. (b): Cumulative distribution functions comparing the three temperature distributions from panel (a). 


the Arecibo Observatory, whose ~ 3.5' beam at 21-cm 
provides the best available single-dish complement to the 
~ 1" VLA beam so that we come closer by orders of mag¬ 
nitude than previous studies to sampling the same gas 
populations in absorption and emission. By obtaining 
extremely high-sensitivity VLA Hi absorption lines, we 
detect signatures of both cold and warm ISM directly in 
absorption and, in combination with emission, directly 
determine the temperatures and column densities of in¬ 
dividual components along the line of sight. We describe 
the observation and data reduction strategies and dis¬ 
cuss details of line fitting and parameter estimation in 
this paper. Important initial results from the first half 
of the survey include: 

1. We achieve median RMS noise in optical depth 
of a T ~ 9 x 10~ 4 per 0.42 km s” 1 channel over 
the 31 sightlines probed by 21-SPONGE so far. 
This allows us to directly probe absorption sig¬ 
nals from both cold and warm Hi. We are able 
to achieve this sensitivity in part by maximizing 
the signal-to-noise in our bandpass calibration so¬ 
lutions, which we accomplish by characterizing pe¬ 
riodic structures within the solutions and combin¬ 
ing calibration observations over time. 

2. Following a careful Gaussian decomposition of all 
31 Hi absorption detections, we estimate T s and 
N(Hi) for 157 Gaussian components. We find a 
maximum spin temperature per individual compo¬ 
nent of T s = 1451 ± 263 K, which is more than 
a factor of 2 higher than the maximum Gaussian 


component spin temperature measured by the Mil¬ 
lennium Arecibo Hi Survey (Heiles & Troland 2003; 
HT03). We are sensitive to lower N(Hi, CNM) due 
to our better sensitivity in optical depth, and de¬ 
tect individual clouds with N(Hi, CNM) as low as 
3± 1 x 10 16 cm -2 . 

3. We compute total CNM column density, total 
WNM column density, total Hi column density 
and CNM fraction along each of the 31 emis¬ 
sion/absorption pair sightlines. We detect simi¬ 
lar total CNM column densities, defined as the to¬ 
tal column density of all Gaussian components per 
sightline with T s < 200 K, as HT03 and Roy et 
al. (2013). We find larger values of total “warm 
gas” column density, defined as the total column 
density of all absorption-detected Gaussian compo¬ 
nents with T s > 200 K and the components only de¬ 
tected in emission, due to the fact that we are sen¬ 
sitive to warmer gas in absorption than HT03. We 
find consistent CNM fractions (/cnm) with those 
found by HT03 and Roy et al. (2013), although the 
Roy et al. (2013) values are all upper limits. 

4. We observe the same trend between CNM fraction 
and harmonic mean LOS spin temperature found 
by the synthetic-observation analysis of a hydro- 
dynamic Galaxy simulation by Kim et al. (2014). 
Although we find lower CNM fractions (median 
0.20) than the simulation predicts (98% LOS with 
0.4 < /cnm < 0.7), the agreement in the trend 
is encouraging given that our calculations rely on 
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Fig. 11.— (a): FT spectrum of a bandpass solution, with the boundaries of the applied step-function filter overlaid as dashed lines, (b): 
Model of the bandpass solution (solid black), computed as the inverse FT of the filter-convolved FT spectrum of the bandpass solution. 
The bandpass solution is overlaid as the dashed red line. 


Gaussian decomposition to estimate CNM frac¬ 
tion and theirs are taken directly from the simula¬ 
tion, suggesting that the Gaussian decomposition 
method is successful. 

5. In a comparison of T s estimation methods, we find 
that one-phase or isothermal assumptions overesti¬ 
mate T s values relative to individual Gaussian fits. 
The T a distribution produced by multi-phase fits 
via Gaussian decomposition better reproduces the¬ 
oretically expected lower limits for individual CNM 
clouds based on detailed heating and cooling ana¬ 
lytical models. 
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APPENDIX 

BANDPASS CALIBRATION 
59 kHz Ripple 

We find that all bandpass solutions obtained by 21-SPONGE contain a sinusoidal “ripple” with a period of 59 kHz 
(or FWHM Av ~ 12 km/s at our velocity resolution) and amplitude of ^0.0015. The period, amplitude and phase of 
this ripple are constant between array configurations and time of year. After investigating the stability of the ripple, 
we determined that it is caused by finite impulse response (FIR) filters in the WIDAR correlator applied prior to 
correlation which are used to shape the bandpass. The constant relative amplitude, phase and period of the ripple can 
all be explained in this context. 

To estimate the quantitative properties of the ripple, we first fit and remove the linear bandpass slope and zoom 
in on the central channels of the bandpass solution (channels 50-200 out of a total 256, shown in Figure mb). The 
linear slope does not affect our ability to detect wide absorption lines, and we are only interested in determining the 
parameters of the ripple and correcting for it. 

We then compute the Fourier transform (FT) spectrum of this flattened, zoomed solution (see Figure II lh ) . Next, 
we multiply this FT spectrum with a simple step function filter to isolate the dominant periodic components. The 
boundaries of the filter are shown in Figure [Tlf i as vertical dashed lines. We then compute the inverse FT of the 
filter-multiplied FT spectrum (i.e., the thick black line in Figure 1 1 lb .) to model the periodic component. The modeled 
signal is displayed by the solid black line in Figure UTb . along with the original slope-subtracted bandpass solution 
(dashed red line). 
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Channels Model residual value 


Fig. 12.— (a): Residuals of the FT model to the bandpass solution displayed in Figure HT1 (i.e.. bandpass solution minus model), (b): 
Histogram of the residuals to check if they are well-represented by a Gaussian distribution, as would be expected in the case of a successful 
model. A Gaussian model is overlaid (solid black line) to illustrate this. 


As a test of the success of the model, we compute the residuals by subtracting the model from the bandpass solution 
(see Figure IT2h). We then compute a histogram of these residuals to observe if they are well-represented by a Gaussian 
distribution, which we expect if the model is reasonable (see Figure [T^b). In this example, a K-S test between the 
residuals and a Gaussian distribution returns p=0.55, indicating that the noise is consistent with being drawn from 
a Gaussian distribution. From this, we can conclude that the model does a reasonably good job of isolating the 
dominant periodic components of the bandpass solution. We compute the RMS noise in these residuals as an indicator 
of the noise level in the bandpass solution that is not due to the periodic components. We record various parameters 
associated with the ripple fit in Table [4] for all bandpass solutions that we tested, including RMS noise in the model 
residuals, ripple amplitude, period and phase. The amplitu de, per iod and p ha se o f the ripple for each observation were 
calculated using a periodogram, following the methods of lHorne fe Baliunasl ([1986D . 

To test if the presence of the ripple is a function of particular antennas, we separated the bandpass solutions into 
three groups of about 9 antennas (i.e antennas 1-9, 10-18, 19-28). After fitting the solutions from these antenna 
subsets, we find that the ripple exists at a similar power as in the full solution in all cases. Therefore, we conclude 
that the effect is not antenna-based. 


Time Averaging 

To minimize noise levels in bandpass (BP) solutions, we combine bandpass observations acquired at different times. 
To compare noise levels in the BP solutions from the original observation and the combined observations, we first 
model and remove the 59 kHz ripple present in all solutions. Figure l~HH shows an example of the modeling process. A 
zoom-in onto a bandpass solution from one observation of 3C147 (33 min integration) is shown in the top left panel, 
with the ripple model overlaid in red. The top right panel displays the residuals following the removal of the ripple 
model, from which the noise in the bandpass solution is computed. The bottom panels are the same, after adding an 
additional observation of 3C147 from 4 days later (for a total of 33+30=63 min integration). 

We expect the noise in the bandpass solution to decline according to the factor of increase in integration time given 
by the ideal radiometer equation, or by 1/■%/^ 2 /^ 1 - For example, from Figure 15, for an increase in integration time of 

ti/ti = 63/33 = 1.9, we expect an improvement in noise by a factor of 1/ s/tvjti = 0.72. The improvement in noise is 
0T0055/0.00073 = 0.75. 

In Figure [lj] we show the RMS noise in the bandpass solution (computed after removing the ripple) as a function of 
for several different observations (symbols). The solid line displays a fit to the data points, which is linear in log 
space with a slope of —0.44 ± 0.02. This is nearly equal to the theoretically expected fit index of —0.5 (e.g. a noise 
decrease by factor of l/i/t^Ti)- This result is highly encouraging, as it demonstrates that VLA bandpass solutions at 
this resolution from different days are stable enough to be combined to increase the sensitivity of observations while 
significantly saving observing time, which is especially useful for very deep observations. 

In this process, we are able to combine observations of different sources. In Figure [M] the green triangles are 
from a trial combining observations of the sources 3C410 («Si.4GHz = 10 Jy) and 3C454.3 (5i.4GHz = llJy) as 
bandpass calibrators. The other three cases are combinations of observations of the same bandpass calibrators, 3C147 
(51.4GHz = 23 Jy) and 3C286 (5i.4GHz = 15 Jy) between several days. These examples represent a common trend in 
our results. 

In later observations, we experimented with including additional 0.5 MHz subbands spaced ±1.5, 3and4.5 MHz from 






























































The 21-SPONGE Survey I 


21 


TABLE 4 

Bandpass Solution Parameters 


Information: 

Target Name 

Obs. Date 

Array Config. 

Fit Parameters: 

RMS noise a 
(per channel) 

Amplitude^ 

Period^ 

(channels) 

Phase c 

(channel) 

3C120D1 

037TT 

B 

0.00049 

0.0014 

30.1 

120 

3C120D2 

03/11 

B 

0.00047 

0.0014 

29.9 

120 

3C286 

05/11 

B 

0.00064 

0.0016 

29.9 

122 

3C225BD1 

05/11 

BnA 

0.00049 

0.0015 

30.0 

121 

3C225BD2 

05/11 

BnA 

0.00049 

0.0015 

30.3 

121 

3C225BD3 

05/11 

BnA 

0.00042 

0.0013 

30.1 

120 

4C12.50 

05/11 

BnA 

0.00048 

0.0014 

30.0 

120 

3C345 

06/11 

BnA-A 

0.00074 

0.0017 

30.1 

120 

3C298 

06/11 

BnA-A 

0.00057 

0.0015 

29.8 

121 

4C32.44D1 

05/11 

BnA-A 

0.00040 

0.0014 

29.9 

122 

4C32.44D2 

06/11 

BnA-A 

0.00063 

0.0015 

29.9 

120 

4C32.44D3 

06/11 

BnA-A 

0.00057 

0.0015 

29.8 

122 

3C48 

08/11 

A 

0.00048 

0.0014 

30.0 

120 

4C16.09 

09/11 

A 

0.00037 

0.0014 

30.0 

120 

4C16.09.2 

09/11 

A 

0.00031 

0.0013 

29.9 

120 

3C133 

09/11 

A-D 

0.00035 

0.0014 

30.1 

120 

PKS0531 

09/11 

A-D 

0.00031 

0.0013 

29.9 

120 

3011 

04/12 

C 

0.00042 

0.0014 

29.9 

120 

3C154D1 

04/12 

C 

0.00057 

0.0017 

30.2 

122 

3C154D2 

04/12 

C 

0.00033 

0.0013 

29.9 

120 

3023 

05/12 

CnB 

0.00041 

0.0015 

30.1 

120 

3038 

05/12 

CnB 

0.00041 

0.0014 

29.9 

122 

3C410 

07/12 

B 

0.00057 

0.0016 

29.6 

123 


a : the RMS noise calculated in the residuals of the bandpass solution after removing the ripple model. 
b : calculated using the periodogram analysis of the solution (i.e. these are the parameters of the dominant component) 
c : starting channel of a new period nearest to the center of the solution. 

the Hi line. We found that combining these subbands does not significantly improve the noise in the bandpass solution. 
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Fig. 13. — An example of combining bandpass solutions. We fit the bandpass solutions (left) for the presence of a 59 kHz ripple, remove 
this model (shown in red) to produce the panels on the right, from which we measure the RMS noise in the solution. The top panels are 
from a single observation, and the bottom panels are the results of adding an observation 4 days later. 
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Fig. 14.— The improvement in RMS noise in a bandpass solution (y-axis) by combining observations in time, thereby increasing the total 
integration time (x-axis). Symbols denote measurements, and the solid line is a linear fit to the data. The index of the fit is —0.44 d= 0.02, 
compared to the theoretical expectation of an index equal to —0.5 (e.g. noise decrease by factor of 1 /shown by the dotted line). 








